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I summarize the current status of the comparison between experiment and the 
predictions of the NRQCD factorization approach to quarkonium decay and pro- 
duction. I also present the results of some recent calculations and theoretical 
developments in the NRQCD factorization approach. 



1. Some Successes of the NRQCD Factorization Approach 

Since its formulation, the NRQCD factorization approach 1 has enjoyed a 
number of successes in predicting inclusive decay rates and production cross 
sections for quarkonium states. 

The initial success of the NRQCD factorization method was in the com- 
putation of infrared-finite predictions for the inclusive decay rates of P-wave 
quarkonium states. 2 A subsequent global fit to the P-wave charmonium 
data 3 has yielded a \ 2 P er degree of freedom of 15.0/10 and values of the 
nonperturbative NRQCD matrix elements that are in good agreement with 
lattice determinations 4 and with estimates based on the NRQCD velocity- 
scaling rules. 1 ' 5 

In the area of quarkonium production, the first success of the NRQCD 
factorization approach was in explaining the Tevatron data for the inclusive 
production of J/ip, Xc, ip(2S), T, and T(2<5) states. Previous calculations 
in the color-singlet model had yielded results that were smaller than the 
data by more than an order of magnitude. In the NRQCD predictions, the 
unknown NRQCD matrix elements were determined from fits to the data. 
Nevertheless, the comparison with the Tevatron data gives a nontrivial 
confirmation of the NRQCD factorization approach because the extracted 
matrix elements satisfy the velocity-scaling rules, and the shape of the data 
as a function of px is consistent with NRQCD factorization, but not with 
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the color-singlet model. 

A recent success of the NRQCD factorization approach is in predicting 
the cross section for 77 — > J/tp + X at LEP. 6 In this case, the predic- 
tion makes use of NRQCD matrix elements extracted from the Tevatron 
data, and hence, provides a test of the predicted universality of the matrix 
elements. The principal theoretical uncertainties arise from uncertainties 
in the renormalization and factorization scales and from estimates of the 
color-octet matrix elements. The data from the Delphi experiment are con- 
sistent with the NRQCD factorization prediction and clearly disfavor the 
prediction of the color-singlet model. 

Another prediction that makes use of the NRQCD matrix elements ex- 
tracted from the Tevatron data is that for quarkonium production in deep- 
inelastic scattering at HERA. 7 Again, the theoretical uncertainties arise 
mainly from uncertainties in the renormalization and factorization scales 
and in the color-octet matrix elements. Data from the HI experiment plot- 
ted as functions of pt or Q 2 favor the prediction of NRQCD factorization 
over that of the color-singlet model. Neither prediction is entirely consis- 
tent with the data as a function of z. It should be noted that the most 
recent calculation 7 disagrees with a number of previous theoretical results, 
and that these disagreements have not yet been resolved fully. 

2. Some Problematic Comparisons with Experiment 

There are some processes for which the comparisons between the experi- 
mental data and the NRQCD factorization predictions are less satisfactory. 

The production of transversely polarized quarkonium at the Tevatron 
is potentially a "smoking gun" for the color-octet production mechanism, 
which is an integral part of the NRQCD factorization formalism. For 
quarkonium production at large px {pt <; 4m c for J/tp) , gluon fragmenta- 
tion into quarkonium via the color-octet mechanism is the dominant pro- 
cess. At large pt, the fragmenting gluon is nearly on mass shell, and, 
so, is transversely polarized. NRQCD predicts 8 that the polarization is 
largely transferred to the J/tp, although it is diluted somewhat by non- 
fragmentation processes, radiative corrections, and feeddown from higher 
quarkonium states. 9,10 ' 11 The data from the CDF experiment show no ev- 
idence for the predicted increase of transverse polarization with increasing 
Pt- However, the error bars are large, and only the highest-pr data point 
is actually inconsistent with the prediction. There are also large theoret- 
ical uncertainties, primarily from corrections of higher order in a and v, 
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including v 2 effects in the polarization transfer. 

Another problematic process is inelastic quarkonium photoproduction 
at HERA. The HERA data 12 ' 13 are just barely compatible with the 
prediction. 14 ^ 15 ' 16 ' 17 The color-octet mechanism leads to a prediction of 
increasing rate with increasing energy fraction z. This is not observed. 
However, theoretical uncertainties arising from corrections of higher order 
in a s and v, uncertainties in m c , and the breakdown of the v expansion near 
z = 1 have not been addressed or have been addressed incompletely. 18 ' 19 ' 20 
It should be noted that the experimental data differential in pt are com- 
patible with color-singlet production alone, 18 even at large pr- 

A recent striking result from the Belle collaboration 21 concerns produc- 
tion of a double cc pair in e + e~ collisions: a(e + e~ — > J/tp cc)/a(e + e~ — > 
J/tpX) = 0.59j:o'.i3±0.12. Pcrturbative QCD plus the color-singlct model 
leads to the prediction 22 cr(e+e~ — » J/i/>cc)/cr(e+e~ — » J/tpX) w 0.1. This 
seems to be a major discrepancy between theory and experiment. 

3. General Difficulties in the Theory 

There are several difficulties that arise when one attempts to make accurate 
theoretical predictions for quarkonium production and decay processes. 

Foremost among these is the fact that the color-octet NRQCD matrix 
elements are poorly determined. For some processes, only linear combina- 
tions of color-octet matrix elements can be determined from the data with 
reasonable accuracy. Different linear combinations arc fixed by different 
processes — a situation that makes it difficult to test the universality of the 
color-octet matrix elements. 

A further difficulty in making accurate theoretical predictions is the fact 
that corrections to the short-distance coefficients of ncxt-to-leading-order 
(NLO) in a s are often large. For example, for the process J/ip — * 777, 
the NLO correction is — 12.62a s /7r, relative to the leading-order coefficient. 
A related problem is that the dependence of the calculated short-distance 
coefficients on the renormalization scale is often large. These issues raise 
doubts about the convergence of the perturbation expansion for the short- 
distance coefficients. For some key processes, such as J/ip production in 
hadronic collisions, the NLO corrections have not yet been computed. For 
certain processes, resummations of large logarithms of p^/M 2 , M|/M 2 , 
and z 2 have been carried out. 23,24 In some instances, resummations of log- 
arithms of 1 — x and x may also be required. 

Another significant source of inaccuracy in theoretical predictions is 
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the existence of large relativistic corrections. For example, for the process 
J/tp — > 777, the order-v 2 correction is — 5.32v 2 . Since v 2 fa 0.3 for char- 
monium and v 2 w 0.1 for bottomonium, one can question whether the w 
expansion converges. Furthermore, for many processes, the order-u 2 cor- 
rections have not yet been computed. 

In the remainder of this paper, I describe some recent progress in con- 
fronting these obstacles to accurate theoretical predictions. 

4. Relativistic Corrections to Gluon Fragmentation into 
S-wave Quarkonium and to 5-wave Quarkonium Decay 

Motivated by the importance of the fragmentation process for S'-wave 
quarkonium production at large pt, J- Lee and I calculated the order- 1> 2 
corrections to that process. Our preliminary results are that the correc- 
tion to the color-singlet process is approximately 74% for the J/tp, and the 
correction to the color-octet process is approximately -54% for the J/tp. 
The correction to the color-singlet process reduces the predicted value of 
the transverse polarization parameter a by about 10% at large px- The 
corrected prediction is still far above the highest-pr CDF data point. The 
correction to the color-octet process directly affects the value of the 3 Si 
color-octet matrix element that is obtained by fitting to the Tevatron data. 

A. Petrelli and I have computed the coefficients of color-singlet opera- 
tors of order w 4 for S-wave quarkonium decays. 25 This is the first decay or 
production calculation at next-to-next-to-leading order in v. The series for 
1 S'o decay into two photons and for the color-singlet part of 1 Sq decay into 
light hadrons are 1 — 1.33u 2 + 1.51v 4 + • • •, the series for 3 Si decay into 

e + e~ is 1 — 1.33v 2 + 1.61u 4 H , and the series for the color-singlet part 

of 3 Si decay into light hadrons is 1 - 5.32w 2 + 7.62v 4 H . All of these 

series seem to be converging in order w 4 , even the one for 3 Si decay into 
light hadrons, which receives a large contribution in order v 2 . It is known 
that color-octet contributions in the first non-trivial order (v 3 and v 4 ) are 
large. 3 ' 25 The hope is that these, too, will receive small corrections in the 
next order in v. 

5. Lattice Computation of Bottomonium Decay Matrix 
Elements 

D.K. Sinclair, S. Kim, and I have recently computed matrix elements for 
bottomonium decay on the lattice, using two dynamical light quarks. 26 
The essential steps in this computation are to measure matrix elements in 
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a lattice simulation, to use perturbation theory (at the cutoff scale mj,) to 
relate the lattice and continuum matrix elements, and to extrapolate to 
three light-quark flavors and to physical light-quark masses. For the ma- 
trix element of the color-singlet 5*-wave operator in the T state, we obtain 
4.10(1)(9)(41) GeV 3 , which compares well with the phenomcnological value 
3.86(14) GeV 3 extracted from the measured rate for T — > e+e~. We also 
obtained values for matrix elements of the color-singlet P-wave operator 
and the color-octet S'-wave operator in the Xb state, which can be used to 
make predictions for the decays of Xb states. The matrix element of the 
order-v 2 color-singlet 5-wave operator in the T state is poorly determined, 
owing to large lattice-to-continuum corrections. The values of all of these 
matrix elements are strongly affected by the inclusion of light dynamical 
quarks in the calculation. 



6. Resummation of QCD Corrections to Quarkonium 
Decay Rates 

In this Section I describe work carried out in collaboration with Y.-Q. 
Chen. 27 

Expressions for the decay rates of the r\ c into light hadrons (two gluons 
in perturbation theory) and two photons are both known through NLO 
in a s and v. Each expression depends on two NRQCD matrix elements. 
However, if we take the ratio of the rates, the dependence on the matrix 
elements cancels: 



/199 13tt 5 

[-6 — 8-;~9 nj 



~ ! h>n JnJ^+0(a 2 s ) + 0(v 3 )^, (1) 



4to 2 

where R Q = 9a 2 s (n)/(8a 2 em ), and /3 = (33 - 2n f )/(6ir). 

One might hope that Eq. (1) would yield an accurate prediction for 
the ratio. Unfortunately, the term in Eq. (1) proportional to a s is ap- 
proximately l.li?o, casting doubt on the convergence of the perturba- 
tion expansion. Ignoring this difficulty and setting fi — 2m c , we obtain 
i? NLO (2m c ) = 2.1 x 10 3 . The experimental value is i? Ex P = (3.3±1.3) x 10 3 . 
The agreement is reasonable, however the \i dependence of the theoreti- 
cal result is strong. For example, using BLM scale setting 28 we obtain 
Mblm ~ 0.52to c . This leads to an NLO term that is 0.5 times the lead- 
ing term, and, so, the convergence appears to be better. However, the 
agreement with experiment is poor: R NLO (hblm) = 9.9 x 10 3 . 
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It is instructive to re- write Eq. (1) in terms of (3q- 
'37 !3tt 2 
~2 



R NLO (») = R 



, . 16 „, [i 2 



(2) 

The terms proportional to a s that do not contain a factor (3 have a value 
2.5i? 0; while the terms proportional to (3 have a value 12i? for [i 2 = 
Am 2 .. The fact that the terms proportional to 0o dominate in the NLO 
correction suggests that we should resum the contributions proportional 
to (a s (3o) n to all orders in a s . To do this, we make use of the method 
of naive non-Abelianization (NNA). 29 That is, we sum the fermion-loop 
vacuum-polarization contributions and then take gluon loops into account 
by replacing the fermion-loop contribution to (3q with the full /3q. 

In order to carry this out, one must compute all cuts through the two 
final-state complete gluon propagators, where each propagator consists of 
a chain of free gluon propagators and vacuum-polarization bubbles. The 
result is 

^ dx_ * 

2irx . 
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(3) 



where x = k 2 /(Am 2 ), y = I 2 /(Am 2 ), k and I are the gluon momenta, 
f(x,y) is the phase-space and heavy-quark factor (/(0,0) = 1), Ir\x) = 
— 21m [— «IT(x)] n is the sum of quark cuts of the complete gluon propagator 
in nth order, Gy = [— ill(0)] n is the sum of gluon cuts of the complete 
gluon propagator in nth order, and H^(x) = (k 2 g^ v — k ll k v )Yl(x) is the 
one- loop vacuum polarization. In Eq. (3), the first, second, and third terms 
correspond to the qq, qg, and gg cuts, respectively. 

The individual terms in Eq. (3) are infrared divergent, but we can write 
R Bub = G 2 + G 2l where 



(4) 



and 



dx 
2nx 



^Lf(x,y)I R {x)^I^\y) 6(1 



(y) 



(5) 
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Gi is infrared finite because of the KLN theorem, and Gi is infrared fi- 
nite because the integrands in its two terms become identical when x or y 
vanishes. 

Let us first consider the properties of a single vacuum-polarization- 
bubblc chain, as manifested in G\. A calculation of G\ yields 

G uw = — r~^~ arctan -j mn 

dx Im [a s ((j,) /3o ( d — lnx + in )]", (6) 



where d is a renormalization-scheme-dependent constant. Now, integra- 
tion of In™ a; down to x = produces factorial growth: J Q dx x m In" x = 
(— l) n n!/(m+ l) n+1 . Therefore, the second term of Gi may contain contri- 
butions that grow as n\. This is an indication that, owing to the failure of 
asymptotic freedom at small x, the perturbation series may diverge. Fac- 
torial growth in the perturbation series is characteristic of the presence of 
a renormalon — a singularity in the Borel transform of the decay rate — that 
signals the importance of nonperturbative effects. 

However, when the imaginary part of Eq. (6) is computed, a seemingly 
miraculous cancellation of the factorial growth occurs. The singularity 
in the Borel transform also vanishes in the imaginary part. In order to 
understand why this happens, let us re-write the second term of G\ as a 
contour integral: 



E 



/ ^I m [a s p Q (-lnx + d + m)} n = y / ^ [a a 0o{- lnx + d)} n , (7) 
Jo f ' „ Jc 2lT 



where the contour runs from —1 to to 1 around the cut in lnx. We 
can deform the contour out of the region of small x into a circle of radius 
unity. Now the entire contour lies in a region in which perturbation theory 
applies. No factorial growth occurs. The series is bounded by a geometric 
series and, hence, is convergent. 

Suppose that, for the second term of G\, we carry out the perturbation 
summation before performing the integration. Then, summing the resulting 
geometric series, we obtain 

1 



L 



1 dx T 
Im 



_l-a s j3 (-lnx + d + in)\' ^ 
This differs from the expression (7) by a small, but nonzero, amount: 
[l/(a s /3o)] exp{-l/[a s (A i )/3o] + d} w e d AQ CD /(a s /3oM 2 ), which has the typ- 
ical form of a nonperturbative contribution. Both orders of operations lead 
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to convergent expressions. However, log a; becomes unbounded at small x, 
and, so, the convergence is not uniform, and different orders of limits can 
produce different results. If we carry out the summation before perform- 
ing the integration, as in the expression (8), there is a Landau pole at in 
the integrand at xq = exp[— l/(a s /3 ) + d}. The residue at this pole is the 
difference between the expressions with different orders of operations. 

This raises a question as to which order of limits, if either, is correct. 
In the unintegrated expression, for x sufficiently small, we are outside the 
radius of convergence of the perturbation series, and the perturbation sum 
is unreliable. In the integrated expression, at any finite order in a s , we can 
deform the x contour to get out of the region of small x. Then, the pertur- 
bation expansion converges, and we can take the n — > oo limit. Hence, this 
latter order of limits is the only one of the two that leads to a perturbation 
expansion that is not obviously unreliable. 

Now let us consider the properties of two vacuum-polarization-bubble 
chains, as manifested in G2 [Eq. (5)]. For x ~ 1, kinematics force y to be 
in the nonperturbative region near 0. But, x ~ 1 implies that \fW is large 
compared with I or with the three-momentum of the heavy quark. Then, in 
NRQCD, the propagator carrying k can be shrunk to a point. The result is 
a one- loop correction to the color-octet operator (D$( 3 Si). We subtract this 
nonperturbative piece from the perturbative calculation and absorb it into 
the quarkonium matrix element of Os( 3 Si). After these subtractions, the 
x and y contours can be deformed out of the nonperturbative region. Note 
that the contour-deformation argument is crucial. Without it, there would 
be contributions from the nonperturbative region in which both x and y 
are small. Such contributions would not have the form of a contribution 
to an NRQCD operator matrix element and would invalidate the NRQCD 
factorization formula. 

Combining NNA resummation with the complete NLO calculation, wc 
obtain i? NNA = (3.01 ± 0.30 ± 0.34) x 10 3 , for a s (2m c ) = 0.247 ± 0.012. 
This is about 50% larger than i? NLO (2m c ) and is in good agreement with 
experiment. The first uncertainty arises from the uncertainty in a s (2m c ). 
The second uncertainty comes from a velocity-scaling estimate of the un- 
known matrix element of Og( 3 Si), which we simply treat as an error of 
order v 3 . In the resummed expression, the uncertainty from the choice of \x 
is reduced to about 11% because the bubble sum is a renormalization-group 
invariant (at leading order in the (3 function), and, hence, only the residual 
non-bubble-sum contribution depends on /1. As is the case with any resum- 
mation, the result depends on the choice of resummation scheme. We can 
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check that our choice is reasonable by comparing with a different resumma- 
tion scheme, namely, the use of the background-field-gauge gluon vacuum 
polarization instead of naive non-Abelianization. This scheme yields a very 
similar result: i? BFG = (3.26 ± 0.31 ± 0.47) x 10 3 . 

One may ask why the BLM method, which is supposed to account for 
vacuum-polarization corrections, yields a different result than NNA resum- 
mation. The BLM method approximately resums the vacuum-polarization 
correction through a choice of scale. A change of scale generates a geomet- 
ric series: a s (//) w as(/x)[l + 2a s (/x)/?o ln(/z///) + ...]. However, in the 
NNA rcsummation, the geometric series in the integrand may not yield a 
geometric series after integration over the gluon virtualities. For example, 
Gi(2m c ) = l + 1.91a s (2m c ) + 2.47a2(2m c ) + 0.97a 3 (2m c )-4.49af(2m c )- 
11.76a® (2m c ) H , which is nothing like a geometric series. 

7. Summary and Discussion 

The NRQCD factorization approach has had significant successes in predict- 
ing the rates for a number of processes. Among these are inclusive P-wave 
quarkonium decays, quarkonium production at the Tevatron, 77 — ► J/ip+X 
at LEP, and quarkonium production in deep-inelastic scattering at HERA. 
Other processes are, so far, more problematic. The predictions of NRQCD 
factorization are unconfirmed for quarkonium polarization at the Tevatron, 
inelastic quarkonium photoproduction at HERA, and double cc production 
at Belle. For the former two processes, owing to large error bars, there is 
no real discrepancy at present. 

More precise theoretical predictions are hampered by uncertainties in 
the NRQCD matrix elements and large corrections in NLO in a s and v, 
which have cast some doubt on the convergence of the series in a s and v. 
Lattice calculations can help to pin down the decay matrix elements. It is 
not yet known how to formulate the calculation of production matrix ele- 
ments on the lattice. Calculations of higher order in v are still lacking for 
many processes, but there is now hope that the v expansion settles down in 
order v 4 . In general, resummation of large corrections to the a s series are 
needed. Standard methods exist for dealing with logarithmic corrections. 
The bubble resummation method is a promising technique for dealing with 
nonlogarithmic corrections associated with vacuum-polarization contribu- 
tions. However, other types of large nonlogarithmic corrections are not yet 
understood. A possible clue to their analysis may lie in the fact that large 
order-a s corrections seem to be correlated with large order-w 2 corrections. 
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It is clear that there are still many interesting and challenging problems 
in heavy-quarkonium physics that remain to be solved. 
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